{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model parameters and mutation effects (EVmutation)\n",
    "\n",
    "## Content\n",
    "\n",
    "This notebook demonstrates how to explore the parameters of undirected graphical models inferred using EVcouplings, and how to use those to quantitatively predict the effects of mutations.\n",
    "\n",
    "## Reference\n",
    "\n",
    "Hopf, T. A., Ingraham, J. B., Poelwijk, F.J., Schärfe, C.P.I., Springer, M., Sander, C., & Marks, D. S. (2017). Mutation effects predicted from sequence co-variation. *Nature Biotechnology* **35**, 128–135 doi:10.1038/nbt.3769\n",
    "\n",
    "## Tutorial\n",
    "\n",
    "### Part 1: Load the model parameters from file\n",
    "\n",
    "This file is generated by the couplings stage of the pipeline (using plmc) and has the extension .model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from evcouplings.couplings import CouplingsModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# load parameters from file to create a pairwise model\n",
    "c = CouplingsModel(\"example/PABP_YEAST.model_params\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 2: Predict mutation effects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from evcouplings.mutate import predict_mutation_table, single_mutant_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read the experimental mutational scanning dataset for PABP by Melamed et al., RNA, 2013\n",
    "data = pd.read_csv(\n",
    "    \"example/PABP_YEAST_Fields2013-singles.csv\", sep=\";\", comment=\"#\"\n",
    ")\n",
    "\n",
    "# predict mutations using our model\n",
    "data_pred = predict_mutation_table(\n",
    "    c, data, \"effect_prediction_epistatic\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# can also add predictions by the corresponding independent model\n",
    "c0 = c.to_independent_model()\n",
    "\n",
    "data_pred = predict_mutation_table(\n",
    "    c0, data_pred, \"effect_prediction_independent\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mutant</th>\n",
       "      <th>linear</th>\n",
       "      <th>log</th>\n",
       "      <th>effect_prediction_epistatic</th>\n",
       "      <th>effect_prediction_independent</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>G126A</td>\n",
       "      <td>0.711743</td>\n",
       "      <td>-0.490571</td>\n",
       "      <td>-2.610615</td>\n",
       "      <td>0.406487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>G126C</td>\n",
       "      <td>0.449027</td>\n",
       "      <td>-1.155127</td>\n",
       "      <td>-5.663638</td>\n",
       "      <td>-0.027602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>G126E</td>\n",
       "      <td>0.588928</td>\n",
       "      <td>-0.763836</td>\n",
       "      <td>-6.611062</td>\n",
       "      <td>-1.827570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>G126D</td>\n",
       "      <td>0.229853</td>\n",
       "      <td>-2.121218</td>\n",
       "      <td>-7.270577</td>\n",
       "      <td>-1.180076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>G126N</td>\n",
       "      <td>0.679435</td>\n",
       "      <td>-0.557593</td>\n",
       "      <td>-5.809167</td>\n",
       "      <td>0.387440</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  mutant    linear       log  effect_prediction_epistatic  \\\n",
       "0  G126A  0.711743 -0.490571                    -2.610615   \n",
       "1  G126C  0.449027 -1.155127                    -5.663638   \n",
       "2  G126E  0.588928 -0.763836                    -6.611062   \n",
       "3  G126D  0.229853 -2.121218                    -7.270577   \n",
       "4  G126N  0.679435 -0.557593                    -5.809167   \n",
       "\n",
       "   effect_prediction_independent  \n",
       "0                       0.406487  \n",
       "1                      -0.027602  \n",
       "2                      -1.827570  \n",
       "3                      -1.180076  \n",
       "4                       0.387440  "
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_pred.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Predicting single-substitution landscape (independent of experiment) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mutant</th>\n",
       "      <th>pos</th>\n",
       "      <th>wt</th>\n",
       "      <th>subs</th>\n",
       "      <th>frequency</th>\n",
       "      <th>column_conservation</th>\n",
       "      <th>effect_prediction_epistatic</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>K123A</td>\n",
       "      <td>123</td>\n",
       "      <td>K</td>\n",
       "      <td>A</td>\n",
       "      <td>0.077201</td>\n",
       "      <td>0.112437</td>\n",
       "      <td>0.796669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>K123C</td>\n",
       "      <td>123</td>\n",
       "      <td>K</td>\n",
       "      <td>C</td>\n",
       "      <td>0.001461</td>\n",
       "      <td>0.112437</td>\n",
       "      <td>-3.337328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>K123D</td>\n",
       "      <td>123</td>\n",
       "      <td>K</td>\n",
       "      <td>D</td>\n",
       "      <td>0.118235</td>\n",
       "      <td>0.112437</td>\n",
       "      <td>-0.316713</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>K123E</td>\n",
       "      <td>123</td>\n",
       "      <td>K</td>\n",
       "      <td>E</td>\n",
       "      <td>0.110503</td>\n",
       "      <td>0.112437</td>\n",
       "      <td>-1.078200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>K123F</td>\n",
       "      <td>123</td>\n",
       "      <td>K</td>\n",
       "      <td>F</td>\n",
       "      <td>0.007791</td>\n",
       "      <td>0.112437</td>\n",
       "      <td>-3.013274</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  mutant  pos wt subs  frequency  column_conservation  \\\n",
       "0  K123A  123  K    A   0.077201             0.112437   \n",
       "1  K123C  123  K    C   0.001461             0.112437   \n",
       "2  K123D  123  K    D   0.118235             0.112437   \n",
       "3  K123E  123  K    E   0.110503             0.112437   \n",
       "4  K123F  123  K    F   0.007791             0.112437   \n",
       "\n",
       "   effect_prediction_epistatic  \n",
       "0                     0.796669  \n",
       "1                    -3.337328  \n",
       "2                    -0.316713  \n",
       "3                    -1.078200  \n",
       "4                    -3.013274  "
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "singles = single_mutant_matrix(\n",
    "    c, output_column=\"effect_prediction_epistatic\"\n",
    ")\n",
    "\n",
    "singles.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Predicting arbitary mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.9225586684769951"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# double mutant L186M, G188A\n",
    "delta_E, delta_E_couplings, delta_E_fields = c.delta_hamiltonian([(186, \"L\", \"M\"), (188, \"G\", \"A\")])\n",
    "delta_E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Short-cuts to single and double mutations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.6052584765675419"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# single mutation matrix\n",
    "c.smm(127, \"E\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6.9225586684769951"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# double mutation matrix (double mutant L186M, G188A)\n",
    "c.dmm(186, 188, \"M\", \"A\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Statistical energy of a sequence (rather than delta to WT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "312.19741128035912"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E, E_couplings, E_fields = c.hamiltonians([c.seq()])[0]\n",
    "E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3: Explore model parameters\n",
    "\n",
    "Please see the documentation of evcouplings.couplings.model.CouplingsModel for all available methods. Note that most of these methods can be accessed using lists of positions/symbols. \n",
    "\n",
    "Examples include:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['K' 'G' 'S' 'G' 'N' 'I' 'F' 'I' 'K' 'N' 'L' 'H' 'P' 'D' 'I' 'D' 'N' 'K'\n",
      " 'A' 'L' 'Y' 'D' 'T' 'F' 'S' 'V' 'F' 'G' 'D' 'I' 'L' 'S' 'S' 'K' 'I' 'A'\n",
      " 'T' 'D' 'E' 'N' 'G' 'K' 'S' 'K' 'G' 'F' 'G' 'F' 'V' 'H' 'F' 'E' 'E' 'E'\n",
      " 'G' 'A' 'A' 'K' 'E' 'A' 'I' 'D' 'A' 'L' 'N' 'G' 'M' 'L' 'L' 'N' 'G' 'Q'\n",
      " 'E' 'I' 'Y' 'V' 'A' 'P' 'H' 'L' 'S' 'R']\n",
      "N\n"
     ]
    }
   ],
   "source": [
    "# full sequence\n",
    "print(c.seq())\n",
    "\n",
    "# symbol for particular position (or list of positions)\n",
    "print(c.seq(127))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140\n",
      " 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158\n",
      " 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176\n",
      " 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194\n",
      " 195 196 197 198 199 200 201 202 203 204]\n"
     ]
    }
   ],
   "source": [
    "# positions in model\n",
    "print(c.index_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['A' 'C' 'D' 'E' 'F' 'G' 'H' 'I' 'K' 'L' 'M' 'N' 'P' 'Q' 'R' 'S' 'T' 'V'\n",
      " 'W' 'Y']\n"
     ]
    }
   ],
   "source": [
    "# alphabet\n",
    "print(c.alphabet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.2060956209897995"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pair coupling parameters\n",
    "c.Jij(127, 172, c.seq(127), c.seq(172))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.30619758"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# field parameters\n",
    "c.hi(127, c.seq(127))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 4: Index mapping for complexes\n",
    "\n",
    "When inferring parameters for a concatenated sequence alignment (e.g. protein complexes or other discontinuous sequence segments), the internal model numbering does not match the numbering of the actual sequence. In this case, the model numbering can be manually remapped such that the model can be indexed using tuples (segment_id, position_in_segment).\n",
    "\n",
    "*(Note this example does not execute here and is intended as a reference of how to use the code only)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from evcouplings.couplings import Segment, SegmentIndexMapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define which segments the concatenated sequence in the model is made of.\n",
    "# Most important parameters are region_start and region_end of Segment\n",
    "s_pare = Segment(\"aa\", \"F7YBW7\", 1, 103, segment_id=\"parE\")\n",
    "s_pard = Segment(\"aa\", \"F7YBW8\", 1, 93, segment_id=\"parD\")\n",
    "\n",
    "# Create index mapper\n",
    "mapper = SegmentIndexMapper(True, 1, s_pare, s_pard)\n",
    "\n",
    "# Update indices in model to complex numbering\n",
    "c_mapped = mapper.patch_model(c, inplace=False)\n",
    "\n",
    "# Now access model using tuples of indices rather than single positions alone\n",
    "print(c.seq((\"parD\", 59)))\n",
    "print(c.smm((\"parD\", 59), \"A\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [python3]",
   "language": "python",
   "name": "Python [python3]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}